! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Copyright (c) 2009, Lawrence Livemore National Security Limited Liability
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are 
! permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this list of 
!    conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice, this list
!    of conditions and the following disclaimer in the documentation and/or other 
!    materials provided with the distribution.
!
! 3. Neither the name of the copyright holder nor the names of its contributors may be 
!    used to endorse or promote products derived from this software without specific prior
!    written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY 
! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 
! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL 
! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 
! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT 
! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
! History
! May 2015 - D. Swales - Modified for COSPv2.0
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MODULE MOD_ICARUS
  USE COSP_KINDS,          ONLY: wp
  USE COSP_PHYS_CONSTANTS, ONLY: amd,amw,avo,grav
  use MOD_COSP_STATS,      ONLY: hist2D
  USE MOD_COSP_CONFIG,     ONLY: R_UNDEF,numISCCPTauBins,numISCCPPresBins,isccp_histTau, &
                                 isccp_histPres
  implicit none
  
  ! Shared Parameters                   
  integer,parameter :: &
       ncolprint = 0 ! Flag for debug printing (set as parameter to increase performance)

  ! Cloud-top height determination
  integer :: &
       isccp_top_height,          & ! Top height adjustment method
       isccp_top_height_direction   ! Direction for finding atmosphere pressure level

  ! Parameters used by icarus
  real(wp),parameter :: &
       tauchk = -1._wp*log(0.9999999_wp), & ! Lower limit on optical depth
       isccp_taumin = 0.3_wp,             & ! Minimum optical depth for joint-hostogram
       pstd = 1013250._wp,                & ! Mean sea-level pressure (Pa)
       isccp_t0 = 296._wp,                & ! Mean surface temperature (K)
       output_missing_value = -1.E+30_wp    ! Missing values

contains
  ! ##########################################################################
  ! ##########################################################################
  SUBROUTINE ICARUS(debug,debugcol,npoints,sunlit,nlev,ncol,pfull,          &
                    phalf,qv,cc,conv,dtau_s,dtau_c,th,thd,frac_out,skt,emsfc_lw,at,&
                    dem_s,dem_c,fq_isccp,totalcldarea, meanptop,meantaucld, &
                    meanalbedocld, meantb,meantbclr,boxtau,boxptop,levmatch)
        
    ! INPUTS 
    INTEGER,intent(in) ::      & !
         npoints,              & ! Number of model points in the horizontal
         nlev,                 & ! Number of model levels in column
         ncol,                 & ! Number of subcolumns
         debug,                & ! Debug flag
         debugcol                ! Debug column flag
    INTEGER,intent(in),dimension(npoints) :: & !
         sunlit                  ! 1 for day points, 0 for night time 
    REAL(WP),intent(in) ::     & !
         emsfc_lw                ! 10.5 micron emissivity of surface (fraction)  
    REAL(WP),intent(in),dimension(npoints) :: & !
         skt                     ! Skin Temperature (K)
    REAL(WP),intent(in),dimension(npoints,ncol,nlev) :: & !
         frac_out                ! Boxes gridbox divided up into subcolumns
    REAL(WP),intent(in),dimension(npoints,nlev) :: & !
         pfull,                & ! Pressure of full model levels (Pascals)
         qv,                   & ! Water vapor specific humidity (kg vapor/ kg air)
         cc,                   & ! Cloud cover in each model level (fraction)
         conv,                 & ! Convective cloud cover in each model
         at,                   & ! Temperature in each model level (K) 
         dem_c,                & ! Emissivity for convective clouds
         dem_s,                & ! Emissivity for stratiform clouds
         dtau_c,               & ! Optical depth for convective clouds
         dtau_s                  ! Optical depth for stratiform clouds 
    REAL(WP),intent(in),dimension(npoints,nlev+1) :: & !
         phalf                   ! Pressure of half model levels (Pascals)!
    integer,intent(in) :: th,thd

    ! OUTPUTS 
    REAL(WP),intent(out),dimension(npoints,7,7) :: & 
         fq_isccp                ! The fraction of the model grid box covered by clouds
    REAL(WP),intent(out),dimension(npoints) :: & 
         totalcldarea,         & ! The fraction of model grid box columns with cloud present
         meanptop,             & ! Mean cloud top pressure (mb) - linear averaging
         meantaucld,           & ! Mean optical thickness
         meanalbedocld,        & ! Mean cloud albedo  
         meantb,               & ! Mean all-sky 10.5 micron brightness temperature
         meantbclr               ! Mean clear-sky 10.5 micron brightness temperature
    REAL(WP),intent(out),dimension(npoints,ncol) :: & 
         boxtau,               & ! Optical thickness in each column
         boxptop                 ! Cloud top pressure (mb) in each column
    INTEGER,intent(out),dimension(npoints,ncol) :: &
         levmatch                ! Used for icarus unit testing only


    ! INTERNAL VARIABLES
    CHARACTER(len=10)                     :: ftn09
    REAL(WP),dimension(npoints,ncol)      :: boxttop
    REAL(WP),dimension(npoints,ncol,nlev) :: dtau,demIN
    INTEGER                               :: j,ilev,ibox
    INTEGER,dimension(nlev,ncol   )       :: acc

    ! PARAMETERS
    character ,parameter, dimension(6) :: cchar=(/' ','-','1','+','I','+'/)
    character(len=1),parameter,dimension(6) :: cchar_realtops=(/ ' ',' ','1','1','I','I'/)
    ! ##########################################################################
    
    call cosp_simulator_optics(npoints,ncol,nlev,frac_out,dem_c,dem_s,demIN)
    call cosp_simulator_optics(npoints,ncol,nlev,frac_out,dtau_c,dtau_s,dtau)

    call ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demIN,skt,emsfc_lw,qv,at,        &
                      pfull,phalf,frac_out,levmatch,boxtau,boxptop,boxttop,meantbclr)

    call ICARUS_COLUMN(npoints,ncol,boxtau,boxptop/100._wp,sunlit,boxttop,&
                       fq_isccp,meanalbedocld,meanptop,meantaucld,totalcldarea,meantb)

    ! ##########################################################################
    ! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM
    ! ##########################################################################
    
    if (debugcol.ne.0) then
       do j=1,npoints,debugcol
          
          ! Produce character output
          do ilev=1,nlev
             acc(ilev,1:ncol)=frac_out(j,1:ncol,ilev)*2
             where(levmatch(j,1:ncol) .eq. ilev) acc(ilev,1:ncol)=acc(ilev,1:ncol)+1
          enddo
          
          write(ftn09,11) j
11        format('ftn09.',i4.4)
          open(9, FILE=ftn09, FORM='FORMATTED')
          
          write(9,'(a1)') ' '
          write(9,'(10i5)') (ilev,ilev=5,nlev,5)
          write(9,'(a1)') ' '
          
          do ibox=1,ncol
             write(9,'(40(a1),1x,40(a1))') &
                  (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev),&
                  (cchar(acc(ilev,ibox)+1),ilev=1,nlev) 
          end do
          close(9)

       enddo       
    end if
    
    return
  end SUBROUTINE ICARUS
  
  ! ############################################################################
  ! ############################################################################
  ! ############################################################################
  SUBROUTINE ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demiN,skt,emsfc_lw,qv,at,        &
                          pfull,phalf,frac_out,levmatch,boxtau,boxptop,boxttop,meantbclr)
    ! Inputs
    INTEGER, intent(in) ::   &
         ncol,               & ! Number of subcolumns
         npoints,            & ! Number of horizontal gridpoints
         nlev                  ! Number of vertical levels
    INTEGER, intent(in), dimension(npoints) :: &
         sunlit                ! 1=day 0=night
    REAL(WP),intent(in) :: &
         emsfc_lw              ! 10.5 micron emissivity of surface (fraction) 
    REAL(WP),intent(in), dimension(npoints) ::  &
         skt                   ! Skin temperature
    REAL(WP),intent(in), dimension(npoints,nlev) ::  &
         at,                 & ! Temperature 
         pfull,              & ! Presure
         qv                    ! Specific humidity
    REAL(WP),intent(in), dimension(npoints,ncol,nlev) :: &
         frac_out,           & ! Subcolumn cloud cover
         dtau,               & ! Subcolumn optical thickness
         demIN                 ! Subcolumn emissivity
    REAL(WP),intent(in), dimension(npoints,nlev+1) :: &
         phalf                 ! Pressure at model half levels

    ! Outputs
    REAL(WP),intent(inout),dimension(npoints) :: &
         meantbclr             ! Mean clear-sky 10.5 micron brightness temperature
    REAL(WP),intent(inout),dimension(npoints,ncol) :: &
         boxtau,             & ! Optical thickness in each column
         boxptop,            & ! Cloud top pressure (mb) in each column
         boxttop               ! Cloud top temperature in each column
    INTEGER, intent(inout),dimension(npoints,ncol)      :: levmatch

    ! Local Variables
    INTEGER :: &
       j,ibox,ilev,k1,k2,icycle
    INTEGER,dimension(npoints) :: &
       nmatch,itrop
    INTEGER,dimension(npoints,nlev-1) :: &
       match
    REAL(WP) :: &
       logp,logp1,logp2,atd
    REAL(WP),dimension(npoints) :: &
       bb,attropmin,attrop,ptrop,atmax,btcmin,transmax,tauir,taumin,fluxtopinit,press,   &
       dpress,atmden,rvh20,rhoave,rh20s,rfrgn,tmpexp,tauwv,wk,trans_layers_above_clrsky, &
       fluxtop_clrsky
    REAL(WP),dimension(npoints,nlev) :: &
       dem_wv
    REAL(WP),dimension(npoints,ncol) :: &
       trans_layers_above,dem,tb,emcld,fluxtop,tau,ptop

    ! ####################################################################################
    ! Compute cloud optical depth for each column by summing up subcolumns
    tau(1:npoints,1:ncol) = 0._wp
    tau(1:npoints,1:ncol) = sum(dtau,dim=3)

    ! Set tropopause values
    if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then 
       ptrop(1:npoints)     = 5000._wp
       attropmin(1:npoints) = 400._wp
       atmax(1:npoints)     = 0._wp
       attrop(1:npoints)    = 120._wp
       itrop(1:npoints)     = 1

       do ilev=1,nlev
          where(pfull(1:npoints,ilev) .lt. 40000._wp .and. &
                pfull(1:npoints,ilev) .gt.  5000._wp .and. &
                at(1:npoints,ilev)    .lt. attropmin(1:npoints))
             ptrop(1:npoints)     = pfull(1:npoints,ilev)
             attropmin(1:npoints) = at(1:npoints,ilev)
             attrop(1:npoints)    = attropmin(1:npoints)
             itrop     = ilev
          endwhere
       enddo

       do ilev=1,nlev
          atmax(1:npoints) = merge(at(1:npoints,ilev),atmax(1:npoints),&
               at(1:npoints,ilev) .gt. atmax(1:npoints) .and. ilev  .ge. itrop(1:npoints))
       enddo
    end if
  
    if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
       ! ############################################################################
       !                        Clear-sky radiance calculation
       !       
       ! Compute water vapor continuum emissivity this treatment follows Schwarkzopf 
       ! and Ramasamy JGR 1999,vol 104, pages 9467-9499. The emissivity is calculated
       ! at a wavenumber of 955 cm-1, or 10.47 microns 
       ! ############################################################################
       do ilev=1,nlev
          press(1:npoints)  = pfull(1:npoints,ilev)*10._wp
          dpress(1:npoints) = (phalf(1:npoints,ilev+1)-phalf(1:npoints,ilev))*10._wp
          atmden(1:npoints) = dpress(1:npoints)/(grav*100._wp)
          rvh20(1:npoints)  = qv(1:npoints,ilev)*amd/amw
          wk(1:npoints)     = rvh20(1:npoints)*avo*atmden/amd
          rhoave(1:npoints) = (press(1:npoints)/pstd)*(isccp_t0/at(1:npoints,ilev))
          rh20s(1:npoints)  = rvh20(1:npoints)*rhoave(1:npoints)
          rfrgn(1:npoints)  = rhoave(1:npoints)-rh20s(1:npoints)
          tmpexp(1:npoints) = exp(-0.02_wp*(at(1:npoints,ilev)-isccp_t0))
          tauwv(1:npoints)  = wk(1:npoints)*1.e-20_wp*((0.0224697_wp*rh20s(1:npoints)* &
                              tmpexp(1:npoints))+(3.41817e-7_wp*rfrgn(1:npoints)))*0.98_wp
          dem_wv(1:npoints,ilev) = 1._wp - exp( -1._wp * tauwv(1:npoints))
       enddo

       fluxtop_clrsky(1:npoints)            = 0._wp
       trans_layers_above_clrsky(1:npoints) = 1._wp
       do ilev=1,nlev
          ! Black body emission at temperature of the layer
          bb(1:npoints) = 1._wp / ( exp(1307.27_wp/at(1:npoints,ilev)) - 1._wp )
          
          ! Increase TOA flux by flux emitted from layer times total transmittance in layers above
          fluxtop_clrsky(1:npoints) = fluxtop_clrsky(1:npoints) + &
               dem_wv(1:npoints,ilev)*bb(1:npoints)*trans_layers_above_clrsky(1:npoints)
          
          ! Update trans_layers_above with transmissivity from this layer for next time around loop
          trans_layers_above_clrsky(1:npoints) = trans_layers_above_clrsky(1:npoints)*&
              (1._wp-dem_wv(1:npoints,ilev))
       enddo

       ! Add in surface emission
       bb(1:npoints) = 1._wp/( exp(1307.27_wp/skt(1:npoints)) - 1._wp )
       fluxtop_clrsky(1:npoints) = fluxtop_clrsky(1:npoints) + &
           emsfc_lw * bb(1:npoints)*trans_layers_above_clrsky(1:npoints)

       ! Clear Sky brightness temperature
       meantbclr(1:npoints) = 1307.27_wp/(log(1._wp+(1._wp/fluxtop_clrsky(1:npoints))))
       
       ! #################################################################################
       !                        All-sky radiance calculation
       ! #################################################################################
       
       fluxtop(1:npoints,1:ncol)            = 0._wp
       trans_layers_above(1:npoints,1:ncol) = 1._wp
       do ilev=1,nlev
          ! Black body emission at temperature of the layer
          bb=1._wp/(exp(1307.27_wp/at(1:npoints,ilev)) - 1._wp)
          
          do ibox=1,ncol
             ! Emissivity
             dem(1:npoints,ibox) = merge(dem_wv(1:npoints,ilev), &
                                         1._wp-(1._wp-demIN(1:npoints,ibox,ilev))*(1._wp-dem_wv(1:npoints,ilev)), &
                                         demIN(1:npoints,ibox,ilev) .eq. 0._wp)

             ! Increase TOA flux emitted from layer
             fluxtop(1:npoints,ibox) = fluxtop(1:npoints,ibox) + dem(1:npoints,ibox)*bb*trans_layers_above(1:npoints,ibox) 
             
             ! Update trans_layer by emitted layer from above
             trans_layers_above(1:npoints,ibox) = trans_layers_above(1:npoints,ibox)*(1._wp-dem(1:npoints,ibox))
          enddo
       enddo

       ! Add in surface emission
       bb(1:npoints)=1._wp/( exp(1307.27_wp/skt(1:npoints)) - 1._wp )
       do ibox=1,ncol
          fluxtop(1:npoints,ibox) = fluxtop(1:npoints,ibox) + emsfc_lw*bb(1:npoints)*trans_layers_above(1:npoints,ibox) 
       end do

       ! All Sky brightness temperature
       boxttop(1:npoints,1:ncol) = 1307.27_wp/(log(1._wp+(1._wp/fluxtop(1:npoints,1:ncol))))

       ! #################################################################################  
       !                            Cloud-Top Temperature
       !
       ! Now that you have the top of atmosphere radiance, account for ISCCP 
       ! procedures to determine cloud top temperature account for partially
       ! transmitting cloud recompute flux ISCCP would see assuming a single layer
       ! cloud. *NOTE* choice here of 2.13, as it is primarily ice clouds which have 
       ! partial emissivity and need the adjustment performed in this section. If it
       ! turns out that the cloud brightness temperature is greater than 260K, then 
       ! the liquid cloud conversion factor of 2.56 is used. *NOTE* that this is 
       ! discussed on pages 85-87 of the ISCCP D level documentation 
       ! (Rossow et al. 1996)
       ! #################################################################################

       ! Compute minimum brightness temperature and optical depth
       btcmin(1:npoints) = 1._wp /  ( exp(1307.27_wp/(attrop(1:npoints)-5._wp)) - 1._wp ) 

       do ibox=1,ncol
          transmax(1:npoints) = (fluxtop(1:npoints,ibox)-btcmin) /(fluxtop_clrsky(1:npoints)-btcmin(1:npoints))
          tauir(1:npoints)    = tau(1:npoints,ibox)/2.13_wp
          taumin(1:npoints)   = -log(max(min(transmax(1:npoints),0.9999999_wp),0.001_wp))
          if (isccp_top_height .eq. 1) then
             do j=1,npoints  
                if (transmax(j) .gt. 0.001_wp .and. transmax(j) .le. 0.9999999_wp) then
                   fluxtopinit(j) = fluxtop(j,ibox)
                   tauir(j) = tau(j,ibox)/2.13_wp
                endif
             enddo
             do icycle=1,2
                do j=1,npoints  
                   if (tau(j,ibox) .gt. (tauchk)) then 
                      if (transmax(j) .gt. 0.001_wp .and. transmax(j) .le. 0.9999999_wp) then
                         emcld(j,ibox) = 1._wp - exp(-1._wp * tauir(j)  )
                         fluxtop(j,ibox) = fluxtopinit(j) - ((1.-emcld(j,ibox))*fluxtop_clrsky(j))
                         fluxtop(j,ibox)=max(1.E-06_wp,(fluxtop(j,ibox)/emcld(j,ibox)))
                         tb(j,ibox)= 1307.27_wp / (log(1._wp + (1._wp/fluxtop(j,ibox))))
                         if (tb(j,ibox) .gt. 260._wp) then
                            tauir(j) = tau(j,ibox) / 2.56_wp
                         end if
                      end if
                   end if
                enddo
            enddo
          endif

          ! Cloud-top temperature
          where(tau(1:npoints,ibox) .gt. tauchk)
             tb(1:npoints,ibox)= 1307.27_wp/ (log(1._wp + (1._wp/fluxtop(1:npoints,ibox))))
             where (isccp_top_height .eq. 1 .and. tauir(1:npoints) .lt. taumin(1:npoints))
                tb(1:npoints,ibox) = attrop(1:npoints) - 5._wp 
                tau(1:npoints,ibox) = 2.13_wp*taumin(1:npoints)
             endwhere
          endwhere
          
          ! Clear-sky brightness temperature
          where(tau(1:npoints,ibox) .le. tauchk) 
             tb(1:npoints,ibox) = meantbclr(1:npoints)
          endwhere
       enddo
    else
       meantbclr(1:npoints) = output_missing_value
    end if

    ! ####################################################################################
    !                           Cloud-Top Pressure
    !
    ! The 2 methods differ according to whether or not you use the physical cloud
    ! top pressure (isccp_top_height = 2) or the radiatively determined cloud top
    ! pressure (isccp_top_height = 1 or 3)
    ! ####################################################################################
    ptop(:,:)=0._wp
    levmatch(:,:)=0
    do ibox=1,ncol
       !segregate according to optical thickness
       if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then  
          
          ! Find level whose temperature most closely matches brightness temperature
          nmatch(1:npoints)=0
          do k1=1,nlev-1
             ilev = merge(nlev-k1,k1,isccp_top_height_direction .eq. 2)        
             do j=1,npoints 
                if (ilev           .ge. itrop(j)     .and. &
                     ((at(j,ilev)  .ge. tb(j,ibox)   .and. &  
                      at(j,ilev+1) .le. tb(j,ibox))  .or.  &
                      (at(j,ilev)  .le. tb(j,ibox)   .and. &
                      at(j,ilev+1) .ge. tb(j,ibox)))) then 
                   nmatch(j)=nmatch(j)+1
                   match(j,nmatch(j))=ilev
                endif
             enddo
          enddo

          do j=1,npoints 
             if (nmatch(j) .ge. 1) then
                k1 = match(j,nmatch(j))
                k2 = k1 + 1
                logp1 = log(pfull(j,k1))
                logp2 = log(pfull(j,k2))
                atd = max(tauchk,abs(at(j,k2) - at(j,k1)))
                logp=logp1+(logp2-logp1)*abs(tb(j,ibox)-at(j,k1))/atd
                ptop(j,ibox) = exp(logp)
                levmatch(j,ibox) = merge(k1,k2,abs(pfull(j,k1)-ptop(j,ibox)) .lt. abs(pfull(j,k2)-ptop(j,ibox)))
             else
                if (tb(j,ibox) .le. attrop(j)) then
                   ptop(j,ibox)=ptrop(j)
                   levmatch(j,ibox)=itrop(j)
                end if
                if (tb(j,ibox) .ge. atmax(j)) then
                   ptop(j,ibox)=pfull(j,nlev)
                   levmatch(j,ibox)=nlev
                end if
             end if
          enddo
       else
          do ilev=1,nlev
             where((ptop(1:npoints,ibox) .eq. 0._wp) .and.(frac_out(1:npoints,ibox,ilev) > 0._wp))
                ptop(1:npoints,ibox)=phalf(1:npoints,ilev)
                levmatch(1:npoints,ibox)=ilev
             endwhere
          end do
       end if
       where(tau(1:npoints,ibox) .le. tauchk)
          ptop(1:npoints,ibox)=0._wp
          levmatch(1:npoints,ibox)=0
       endwhere
    enddo

    ! ####################################################################################
    !                Compute subcolumn pressure and optical depth
    ! ####################################################################################
    boxtau(1:npoints,1:ncol)  = output_missing_value
    boxptop(1:npoints,1:ncol) = output_missing_value
    do ibox=1,ncol
       do j=1,npoints 
          if (tau(j,ibox) .gt. (tauchk) .and. ptop(j,ibox) .gt. 0._wp) then
             if (sunlit(j).eq.1 .or. isccp_top_height .eq. 3) then
                boxtau(j,ibox) = tau(j,ibox)
                boxptop(j,ibox) = ptop(j,ibox)!/100._wp
             endif
          endif
       enddo
    enddo

  end SUBROUTINE ICARUS_SUBCOLUMN

  ! ######################################################################################
  ! SUBROUTINE icarus_column
  ! ######################################################################################
  SUBROUTINE ICARUS_column(npoints,ncol,boxtau,boxptop,sunlit,boxttop,fq_isccp,     &
                           meanalbedocld,meanptop,meantaucld,totalcldarea,meantb)
    ! Inputs
    INTEGER, intent(in) :: &
         ncol,    & ! Number of subcolumns
         npoints    ! Number of horizontal gridpoints
    INTEGER, intent(in),dimension(npoints) :: &
         sunlit     ! day=1 night=0 
    REAL(WP),intent(in),dimension(npoints,ncol) ::  &
         boxttop,  & ! Subcolumn top temperature
         boxptop,  & ! Subcolumn cloud top pressure
         boxtau      ! Subcolumn optical depth

    ! Outputs
    REAL(WP),intent(inout),dimension(npoints) :: &
         meanalbedocld, & ! Gridmean cloud albedo
         meanptop,      & ! Gridmean cloud top pressure (mb) - linear averaging
         meantaucld,    & ! Gridmean optical thickness
         totalcldarea,  & ! The fraction of model grid box columns with cloud present
         meantb           ! Gridmean all-sky 10.5 micron brightness temperature 
    REAL(WP),intent(inout),dimension(npoints,7,7) :: &
         fq_isccp         ! The fraction of the model grid box covered by clouds

    ! Local Variables
    INTEGER :: j,ilev,ilev2
    REAL(WP),dimension(npoints,ncol) :: albedocld
    LOGICAL, dimension(npoints,ncol) :: box_cloudy

    ! Variables for new joint-histogram implementation
    logical,dimension(ncol) :: box_cloudy2

    ! ####################################################################################
    !                           Brightness Temperature
    ! ####################################################################################
    if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
       meantb(1:npoints)=sum(boxttop,2)/ncol
    else
       meantb(1:npoints) = output_missing_value
    endif

    ! ####################################################################################
    !                 Determines ISCCP cloud type frequencies
    !
    ! Now that boxptop and boxtau have been determined, determine amount of each of the 
    ! 49 ISCCP cloud types. Also compute grid box mean cloud top pressure and 
    ! optical thickness.  The mean cloud top pressure and optical thickness are 
    ! averages over the cloudy area only. The mean cloud top pressure is a linear
    ! average of the cloud top pressures. The mean cloud optical thickness is 
    ! computed by converting optical thickness to an albedo, averaging in albedo 
    ! units, then converting the average albedo back to a mean optical thickness.  
    ! ####################################################################################

    ! Initialize
    albedocld(1:npoints,1:ncol)  = 0._wp
    box_cloudy(1:npoints,1:ncol) = .false.
    
    ! Reset frequencies
    !fq_isccp = spread(spread(merge(0._wp,output_missing_value,sunlit .eq. 1 .or. isccp_top_height .eq. 3),2,7),2,7)
    do ilev=1,7
       do ilev2=1,7
          do j=1,npoints ! 
             if (sunlit(j).eq.1 .or. isccp_top_height .eq. 3) then 
                fq_isccp(j,ilev,ilev2)= 0._wp
	     else 
                fq_isccp(j,ilev,ilev2)= output_missing_value
             end if
          enddo
       enddo
    enddo

    
    ! Reset variables need for averaging cloud properties
    where(sunlit .eq. 1 .or. isccp_top_height .eq. 3)
       totalcldarea(1:npoints)  = 0._wp
       meanalbedocld(1:npoints) = 0._wp
       meanptop(1:npoints)      = 0._wp
       meantaucld(1:npoints)    = 0._wp
    elsewhere
       totalcldarea(1:npoints)  = output_missing_value
       meanalbedocld(1:npoints) = output_missing_value
       meanptop(1:npoints)      = output_missing_value
       meantaucld(1:npoints)    = output_missing_value
    endwhere
    
    ! Compute column quantities and joint-histogram
    do j=1,npoints 
       ! Subcolumns that are cloudy(true) and not(false)
       box_cloudy2(1:ncol) = merge(.true.,.false.,boxtau(j,1:ncol) .gt. tauchk .and. boxptop(j,1:ncol) .gt. 0._wp)

       ! Compute joint histogram and column quantities for points that are sunlit and cloudy
       if (sunlit(j) .eq.1 .or. isccp_top_height .eq. 3) then 
          ! Joint-histogram
          call hist2D(boxtau(j,1:ncol),boxptop(j,1:ncol),ncol,isccp_histTau,numISCCPTauBins, &
               isccp_histPres,numISCCPPresBins,fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins))
          fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins) = &
               fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins)/ncol
          
          ! Column cloud area
          totalcldarea(j) = real(count(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin),wp)/ncol
             
          ! Subcolumn cloud albedo
          !albedocld(j,1:ncol) = merge((boxtau(j,1:ncol)**0.895_wp)/((boxtau(j,1:ncol)**0.895_wp)+6.82_wp),&
          !     0._wp,box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin)
          where(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin)
             albedocld(j,1:ncol) = (boxtau(j,1:ncol)**0.895_wp)/((boxtau(j,1:ncol)**0.895_wp)+6.82_wp)
          elsewhere
             albedocld(j,1:ncol) = 0._wp
          endwhere
          
          ! Column cloud albedo
          meanalbedocld(j) = sum(albedocld(j,1:ncol))/ncol
          
          ! Column cloud top pressure
          meanptop(j) = sum(boxptop(j,1:ncol),box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin)/ncol
       endif
    enddo
    
    ! Compute mean cloud properties. Set to mssing value in the event that totalcldarea=0
    where(totalcldarea(1:npoints) .gt. 0._wp)
       meanptop(1:npoints)      = 100._wp*meanptop(1:npoints)/totalcldarea(1:npoints)
       meanalbedocld(1:npoints) = meanalbedocld(1:npoints)/totalcldarea(1:npoints)
       meantaucld(1:npoints)    = (6.82_wp/((1._wp/meanalbedocld(1:npoints))-1._wp))**(1._wp/0.895_wp)
    elsewhere
       meanptop(1:nPoints)      = output_missing_value
       meanalbedocld(1:nPoints) = output_missing_value
       meantaucld(1:nPoints)    = output_missing_value
    endwhere
    !meanptop(1:npoints)      = merge(100._wp*meanptop(1:npoints)/totalcldarea(1:npoints),&
    !                                 output_missing_value,totalcldarea(1:npoints) .gt. 0)
    !meanalbedocld(1:npoints) = merge(meanalbedocld(1:npoints)/totalcldarea(1:npoints), &
    !                                 output_missing_value,totalcldarea(1:npoints) .gt. 0)
    !meantaucld(1:npoints)    = merge((6.82_wp/((1._wp/meanalbedocld(1:npoints))-1.))**(1._wp/0.895_wp), &
    !                                 output_missing_value,totalcldarea(1:npoints) .gt. 0)

    ! Represent in percent
    where(totalcldarea .ne. output_missing_value) totalcldarea = totalcldarea*100._wp
    where(fq_isccp     .ne. output_missing_value) fq_isccp     = fq_isccp*100._wp
    
    
  end SUBROUTINE ICARUS_column
  
  subroutine cosp_simulator_optics(dim1,dim2,dim3,flag,varIN1,varIN2,varOUT)
    ! INPUTS
    integer,intent(in) :: &
         dim1,   & ! Dimension 1 extent (Horizontal)
         dim2,   & ! Dimension 2 extent (Subcolumn)
         dim3      ! Dimension 3 extent (Vertical)
    real(wp),intent(in),dimension(dim1,dim2,dim3) :: &
         flag      ! Logical to determine the of merge var1IN and var2IN
    real(wp),intent(in),dimension(dim1,     dim3) :: &
         varIN1, & ! Input field 1
         varIN2    ! Input field 2
    ! OUTPUTS
    real(wp),intent(out),dimension(dim1,dim2,dim3) :: &
         varOUT    ! Merged output field
    ! LOCAL VARIABLES
    integer :: j
    
    varOUT(1:dim1,1:dim2,1:dim3) = 0._wp
    do j=1,dim2
       where(flag(:,j,:) .eq. 1)
          varOUT(:,j,:) = varIN2
       endwhere
       where(flag(:,j,:) .eq. 2)
          varOUT(:,j,:) = varIN1
       endwhere
    enddo
  end subroutine cosp_simulator_optics
end module MOD_ICARUS
